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Abstract 

Effects of non-thermal high-energy electrons on Langmuir wave-particle interaction are investi- 
gated by an initial value approach. A Vlasov-Poisson simulation is employed which is based on the 
splitting scheme by Cheng and Knorr [Cheng, C.Z. and G. Knorr, 1976: J. Comput. Phys. 22, 
330-351.]. The kappa distribution function is taken as an example of non-thermal electrons. The 
modification is manifested as an increase in the Landau damping rate and a decrease in the real 
frequency for a long wavelength limit. A part of the analyses by the modified plasma dispersion 
function [Summers, D. and R.M.Thorne, 1991: Phys. Fluids, B 3, 1835-1847.] is reproduced for 
K = 2,3 and 6. The dispersion relation from the initial value simulation and the plasma dispersion 
function compare favorably. 

PACS numbers: 52.35.Fp, 52.35. Sb, 52.65.Ff 
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I. INTRODUCTION 



Space plasma is far from being in a thermal equilibrium. Suprathermal electrons are 
often observed in space plasmas (Vasyliun 1968). The kappa distribution function (Leubner 
2004) is one of the good examples of non-Maxwellian distribution functions. In one limit 
(Tsallis 1988) the kappa distribution function evolves toward Maxwellian. 

To investigate wave-particle interaction, or Landau damping (Landau 1946; Jackson 1960) 
one needs to incorporate the velocity space dynamics by the Vlasov equation. One of the 
first pieces of work which made the numerical simulation of Vlasov equation available is 
the splitting scheme by Cheng and Knorr (Cheng 1976), which is based on the method of 
characteristics. This method has become a standard method for the Vlasov type simulation. 
The method applies as long as the system is dissipation-less, in other words, if the phase 
volume of the system conserves. The Vlasov-type simulation in lower dimensional cases 
has advantages over Particle-in-Cell (PIC) simulation, since the Vlasov simulation does 
not accompany statistical errors. The Vlasov simulation in lower dimension is suitable for 
investigating subtle effects such as a slight deviation of the equilibrium distribution function 
from Maxwellian. 

One of our long term goals is to investigate the dynamics of Langmuir solitons (Zakharov 
1972) which can then evolve into Langmuir turbulence (Wang 1994, 1995, and 1996). By 
the splitting scheme (Cheng 1976), the electrostatic Vlasov simulation has revealed heating 
of plasmas by Langmuir solitons (Li 1995). 

The purpose of this paper is two-fold. One purpose is to recapitulate the method of 
(Cheng 1976) accurately for the further advanced study of Langmuir solitons. In this work, 
the splitting scheme is revisited and the results of (Cheng 1976) are reproduced. The other 
purpose is, as an initial exercise, to capture the effect of non-Maxwellian distribution function 
on Landau damping. 

This paper is organized as follows. In SecJTll the basic computation model is described. 
In Sec JIIIl we then start from verifying the simulation results with free streaming case whose 
analytical solutions are known. The benchmark linear and nonlinear numerical simulation 
results are discussed in Sec JIVI The effects of non-Maxwellian kappa distribution function are 
discussed in Sec. |Vl Direct comparison of simulation results with modified plasma dispersion 
function is discussed in Sec. I VII We summarize this work in Sec. IVIII 
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II. MODEL EQUATIONS AND NUMERICAL METHODS 



In this section, the model equation of the Vlasov-Poisson simulation is described. For the 
transparency of the work, we recapitulate Cheng (1976) as precise as possible including the 
notations. A one- dimensional Vlasov-Poisson system in the MKS unit is given by 

sA + "&A + ^5;/, = o. (1) 

and 

dE e , . 

— = —{m- raj , (2) 

ox Eo 

where the densities nj = nj (x, t) of each species (subscripts j = i for the ions and j = e for 
the electrons) are given by the distribution functions fj = fj{x,v,t) 



nj{x,t) = / fj{x,v,t)dv. (3) 

J — oo 

Here rrij and qj are the mass and the charge of the species. The electric field is given by E 
and the vacuum permittivity is given by Eq. The coordinates in configuration and velocity 
spaces are given by x and v, respectively. 

1/2 

The equations are further normalized by the Debye length Ae = {EoTe/norrie) , the 

1 /2 ~ 

plasma frequency ooe = {riQe^ / Eorrie) , and the electrostatic field is normalized hj E = 
{eX(.E/Tg) which is equivalent to having the electrostatic energy being comparable to the 
electron thermal energy. Here, e is the unit charge and the electron temperature is given by 
Tf.. By employing the bars denoting the normalized values, we have x = x/Ag, t = oOet, v = 
\/XeUJe, and / = {XeUJe/no)f , where uq is the equilibrium electron and ion densities. Note 
that the thermal velocity of the electrons are given by Ve = (Te/rrieY^'^ = XeCOe- 

After the normalization, we obtain a Vlasov-Poisson system in the Ae (space) and cUe 
(time) scales. 



|/ + V, t) - Epix, V, t) = 0, (4) 



and 



f{x,v,i)dv, (5) 



dx 

where the ion density is taken to be uniform [signified by unity in Eq.(l5])] and the distribution 
function / is for the electrons. Equations (4) and (5) correspond to Eqs.(la) and (lb) of 
Cheng (1976). Hereafter we drop the bars. 
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The splitting scheme (Cheng 1976) is based on the method of characteristics which is 
equivalent to the lowest order symplectic integrator (Ruth 1983). We evolve the distribution 
function by tracing the characteristic curves in the phase space. The method takes three 
steps which is given by 

rix,v) = rix-vAt/2,v), (6) 
r{x,v) = r{x,v + E{x)At), (7) 

and finally, 

r+'{x,v) = rix-vAt/2,v), (8) 

whose kinetic energy part and the potential energy part are time advanced alternatively (the 
lowest order method corresponds to the well-known leap frog method employed frequently 
in PIC simulation). The superscript n stands for the time step. 

Note that the splitting scheme is not a finite difference scheme. In the splitting scheme, if 
the reference points along the characteristic curves "x — or + E{x)At" are exactly 
on the mesh points, the method is quite trivial. However, in general, the points of references 
are located in between mesh points of the x and v space. We thus need an interpolation 
technique to realize Eqs.(l6]),(I7j), and (JHl). As in Cheng (1976), we use Fourier interpolation 
in the configuration space and linear interpolation in the velocity space (Watanabe 2005). 
The computational mesh we employed is exactly that of Fig.l in Cheng (1976). Note that 
we do not have mesh points on the v = axis. 

To make the simulation self-consistent we need to solve Poisson equation, Eq.(l5]). The 
Poisson equation is solved by the Fourier transform since we adopt a periodic configuration 
in X in this paper. All the calculations in this paper employ periodic boundary conditions. 
The electric field is solved directly as in Eq.(|5]) without calculating the electrostatic potential. 



III. FREE STREAMING CASE WITH E = 



To begin with the numerical simulation, we start by verifying the solution of free stream- 
ing case [the Van Kampen mode, (Van Kampen 1955)] whose analytical solutions are known. 
Setting = in Eq.(|4]), we obtain 

-f{x,v,t) + —f{x,v,t) = 0. (9) 
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Following Nicholson 1992, for example, the analytical solution can be given by setting 

f{x,v,t)r^exp{—ikvt). (10) 

For example, if we take an initial condition 

f{x,v,0) = Acos{kx)e-^^/^ (11) 

we obtain the analytical solution 

/(x, v,t) = A cos {kx)e-''^^^ cos (kvt). (12) 

Shown in Fig. 1(a) is the solution of Vlasov equation in a free-streaming case. No force 
is acting on the electrons. The thick solid line is the initial condition of the distribution 
function at a fixed point x = 0. The time advanced distribution function by the numerical 
simulation is given by the black dots at t = 12.5 (time is normalized by the inverse of plasma 
frequency, cj^T^), which matches with the analytical solution given by the dash-dotted curves. 
In the simulation, parameters employed are the maximum cut-off velocity v^ax = lOt'e, and 
< X < 47r. For the mesh points, 32 and 256 are taken in the x and v direction, respectively 
(although we did calculate, the regions |f | > 5.0 are not shown in the figure). Shown in 
Fig. 1(b) is the solution of Vlasov equation in the free-streaming case but with an initial 
condition given by 

f{x, v,0) = [l + A cos {kx)]e-'''/^. (13) 

The time advanced distribution function at t = 6.25 is given by the solid curve at x = vr and 
by the dash-dotted curve at x = 37r. Note that, as demonstrated by the two solutions at 
X = vr and x = 37r, the evolution of the distribution function exhibits point reflection across 
X = 27r. The distribution function streams in positive direction in f > while streams in 
negative direction in f < 0, as a reminder. In both Fig.l(a) and Fig. 1(b), A = 0.5 and 
k = 0.5 are taken. The calculation discussed in this section validates the interpolation 
scheme employed for the Vlasov equation. 

IV. BENCHMARK OF LINEAR AND NONLINEAR SIMULATION RESULTS 

In this section, linear and nonlinear simulation results are compared and benchmarked 
with those of Cheng (1976). For both the linear and the nonlinear simulation, we take an 
initial condition of the form of Eq. ( I13p . 
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In the linear simulation, parameters employed are exactly those of Cheng (1976): the 
maximum cut-off velocity Vmax = 4.0we and < x < An, and the mesh points are 8 and 32 
for X and v, respectively. In Fig.2, the Landau damping phase is shown in terms of electric 
field strength \E\. As in Cheng (1976), recurrence effect takes place after t = 42. Note that 
only the Fourier component k = 0.5 is kept and the other modes are filtered out in the linear 
calculation. A mesh point x = tt is chosen for a diagnostic point in Fig.2. 

The figure corresponds to Fig.3 of Cheng (1976) except that A = 0.01 instead oi A = 0.5 
is taken. The measured frequency and the damping rates are u = 1.41 and 7 = —0.155, 
respectively. 

The nonlinear simulation is shown in Fig.3. In the simulation of Fig.3, A = 0.5 and 
k = 0.5 are taken. Figure 3(a) shows the time evolution of electric field at a fixed point 
X = IT. Instead of monotonic decrease the saturation of the amplitude can be seen after 
t = 20. Figure 3(b) shows the distribution function at t = 75 at a fixed point x = ir. 
In Fig(b), we can see a local flattening of the distribution function in the vicinity of the 
phase velocity. The phase velocity of the Langmuir wave estimated by the linear theory is 
u/k = 2.82. 

Note that in the nonlinear simulation, as suggested in Cheng (1976), the frequencies of all 
the higher modes come into play at the later stage. As a result, resonance occurs at multiple 
locations in the velocity space and thus microscopic structures are generated [manifested as 
wrinkles in Cheng (1976)] whose size can be comparable to the mesh size in the velocity 
space. To resolve all the resonance, one needs to employ an extremely high resolution in the 
velocity space. 



V. EFFECTS OF KAPPA DISTRIBUTION FUNCTIONS 



In this section, we investigate the effects of high energy electrons by employing kappa 
distribution functions (Leubner 2004) instead of a Maxwellian (for the initial condition). A 
kappa distribution function we employed is given by 

,2 ■ 



1 + - 
2k 



-K-l 



fv (v) OC 

Note that Maxwellian and kappa distribution functions are related 

lim /„(w) OC exp(— 1;^/2). 



(14) 



(15) 
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The spatial distribution is given in the form of Eq. (fT3|) for the initial condition, thus 
/(x, V, 0) ~ [1 + A cos {kx)]fy{v) is given. We have normalized the kappa distribution func- 
tion to satisfy J^^dvf^{v) = 1, so that the effective number of electrons will be the same 
as in the Maxwellian case. Some of the notable features of the kappa distribution function 
are shown in Fig.4. In Fig. 4(a), Maxwellian (black) and kappa distributions with n = 1 
(red) are compared. One can see large population of high-energy tail in the k = 1 case. The 
functions are plotted in the logarithmic scales for k, = 1 (red) , k = 2 (green), k = 5 (blue), 
and Maxwellian (black) cases in Fig.4(b). 

Figure 5 shows the linear damping with cases when k, = 1 (red), k, = 2 (green), and 
Maxwellian (black) are taken as initial distribution functions. Parameters employed are 
Vmax = lOfg and < x < 47r with 32 and 256 mesh points in x and v. As in Fig. 2, we have 
taken A = 0.01 and k = 0.5. 

With the kappa distribution function the Landau damping rate increases and the real 
frequency decreases (smaller k has larger effects). The variation of the linear damping rate 
7 and the real frequencies u are summarized in Fig.6(a) and Fig.6(b) as a function of n. 
The value k varies from 1.0 to 10.0. The values plotted in both Fig. 6(a) and Fig.6(b) are 
normalized by that of Maxwellian (thus in the figures, 'y/\'~fMax\ = ~1 and u/uMax = 1 for 
a Maxwellian initial condition). 

From the linear theory (Landau 1946; Jackson 1960) the Landau damping rate is pro- 
portional to the slope of the distribution function at the phase velocity of the wave [see, for 
example, Nicholson (1992)]. Employing the measured real frequencies u of Fig.6(b) (and 
k = 0.5 for the wave number), the slope of the distribution function dyfv{v) at the phase 
velocity u/k is estimated and shown in Fig. 7. Smaller k cases have more negative values 
(and thus larger damping rate) which supports the nature of Fig.6(a). For small values 
of k, these latter trends (damping rate increasing and the real frequency decreasing) are 
consistent with the analytical work (Chateau 1991; Summers 1991; Thorne 1991). In the 
next section, we conduct a direct comparison of the Vlasov simulation with the roots of the 
plasma dispersion relation, by taking exactly the same distribution functions employed in 
Summers (1991) and Thorne (1991). 
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VI. DIRECT COMPARISON WITH MODIFIED PLASMA DISPERSION FUNC- 
TION 

By an analogy of plasma dispersion function (Fried 1961) for Maxwellian, Summers and 
Thorne (Summers 1991; Thorne 1991) have extended their work to kappa distribution func- 
tion. We compare the damping rate and real frequencies in our simulation with their theo- 
retical work for different values of n and different wave-numbers k. 

To see the match between the two, in this section we take exactly the same distribution 
function employed in Summers and Thorne (Summers 1991; Thorne 1991). We have taken 
initial distribution function in the form (see Appendix) 



Note the difference in the exponent part ("— k" dependence instead of "— k — 1") between 
Eq.(IT6D and Eq.([14]). 

In Fig. 8, we plot the simulation results and the numerical roots of the dispersion relation 
for K = 2, 3, and 6 [reproduced from the modified dispersion function of Summers (1991) 
and Thorne (1991)]. Figure 8(a) is for the damping rates 7 versus the wave-numbers k. 
Figure 8(b) is for the real frequencies u versus k. The roots from the (modified) dispersion 
relation are plotted as dash-dotted curves. The black, red, and green dash-dotted curves are 
for K = 2, 3, and 6, respectively. 

Those obtained from linear numerical simulation are plotted as black circles. The damp- 
ing rate and the real frequencies are obtained from oscillation signal of the electric field at 
a fixed point x = vr. In this section we did the survey only up to A; = 2.0. When the magni- 
tudes of linear damping rate and the real frequencies become comparable, measurement of 7 
and u becomes troublesome. The dispersion relation from the initial value simulation com- 
pare favorably with that from the plasma dispersion function. With the smaller k values, 
the real frequency decreases. Note that, however, with the smaller k values, the absolute 
value of the damping rates can also decrease for larger values of /c, contrary to what we have 
obtained in the previous section. The damping rates are sensitive function of the local value 
dyf where the resonant phase velocities "w/A;" are located at, and can vary depending on 
the wave-numbers k. 

One of the advantages of the initial value approach is its application to nonlinear simula- 




(16) 
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tion. Our preliminary nonlinear simulation results employing a kappa distribution function 
as an initial condition are presented below. Figure 9 shows electron distribution functions 
suggesting long time evolution up to t = 1500. The distribution functions are given at a 
fixed point x = n. Figure 9(a) is for a Maxwellian and Fig.9(b) is for a kappa distribution 
function {n = 2). The dash-dotted curves are for t = and the solid curves are for t = 1500. 
The integration of the distribution functions over the velocity space is the same for the two 
cases. The integration of the ion density over the configuration space is kept the same with 
the electron density (total numbers of ions and electrons in the system are the same). In 
Fig.9, a relatively large cut-off velocity v/vthe = 12.0 is taken with a high resolution (1024 
mesh points are taken in the velocity space). The k = 2 distribution function have larger 
population at the high energy tail. In the simulation of Fig.9, A = 0.5 and k = 0.5 are taken. 
Figure 9(c) shows local expansion of Fig. 9 (a) and Fig. 9(b) near the resonant phase veloci- 
ties. In the figure, the perpendicular lines are suggesting the phase velocities. Note that the 
dash-dotted black line tends to solid black line (frequency down-shift for Maxwellian) while 
the dash-dotted red line tends to solid red line (frequency up-shift for kappa function). In 
a very long time scale, normal mode frequency changes and the distribution functions can 
possibly evolve toward a similar equilibrium state. 

VII. SUMMARY 

In this work, the splitting scheme is revisited and the simulation results are compared with 
Cheng (1976). Based on the validation of the code, as an initial exercise, we have discussed 
the effect of non-Maxwellian distribution function by employing the kappa-distribution func- 
tion. 

The slope [the absolute value of dyf{v)] of the distribution function at the phase velocity 
is estimated which supports the nature of increasing damping rate at smaller n values. The 
simulation results compare favorably with the analyses based on modified plasma dispersion 
function (Summers 1991; Thorne 1991). The specific calculations we have demonstrated are 
for K = 2, 3, and 6, with the wave-numbers k = 0.5, 1.0 and 2.0. Our preliminary nonlinear 
simulation employing the kappa distribution function is presented. 

A part of this work is supported by National Cheng Kung University Top University 
Project and a part by National Science Council of Taiwan, NSC 100-2112-M-006-021-MY3. 
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Appendix A: A brief review of modified plasma dispersion function 

We review the modified plasma dispersion function (the "Z* function" ) for kappa distri- 
bution functions (Summers 1991). In this appendix, we invert the notation (see SecHII); the 
values with bars are normalized ones, and all other values are those before normalization. 
Neglecting the ion contribution and requiring the longitudinal component of the dielectric 
tensor to be zero (Fried 1961), we obtain the Langmuir wave dispersion relation 

r dv^^ = (Ai) 

Here, we have assumed a planer wave of the form exp {ikx — iut) for all the perturbed 
quantities, where i is the imaginary unit. Equation flAip is exactly what we have in our 
Vlasov system employed for the numerical simulation. For a one dimensional Maxwellian 

/.(/) = ^exp(-t;V^;D' (A2) 
by introducing the plasma dispersion function (the "Z function". Fried 1961) 

1 f°° p~^^ 

Z{i) = — -ds, (A3) 



we arrive at the well-known Langmuir wave dispersion relation 

2a;? 



1 + 



k'^v'f 



0, (A4) 



where v"^ = 2Te/m.e = 2Wg. 

In contrast. Summers (1991) employs a one dimensional kappa distribution function 



^" ^""^ ~ V^v, V^^TiK - 1/2)^3/2 + 2«: - a) • ^^^^ 

Here, = e~'^t'^~^ds stands for Gamma function. Substituting Eq. flASD into Eq. flAip 
and by introducing a modified plasma dispersion function, 
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we arrive at the dispersion relation for kappa distribution function 



1 + 



2ul K 



=0 



(AT) 



where 



UJ 



\/2k - 3kve' 

Note that the exponent "— k — 1" in Eq.f lA6P arises from taking a derivative on the right 
hand of Eq. flASI) . A normahzed dispersion relation is given by 



where k = kX^.. 

Following Thorne (1991), the root finding algorithm employed in Fig.8 is stated. We first 
let ^ = X + iy, where x and y are real numbers. We then fix the y value and solve the 
imaginary part of the dispersion relation Eg. (1A8I) . Im [^Z*{^)] = 0, to obtain x. When both 
X and y are given, we solve the real part of Eq. (lASP for k by 



+ 1 



1 , ^-^h. 







(A8) 



2ac k 




Finally, the real frequencies and the damping rates are given by a; = kx^y {2k — 3)/k and 



7 = —ky^y (2k — 3)/k. As in SecJTTl the time scale is normalized by ^. The black, red, 
and green curves in Fig.8 are obtained by this latter algorithm. 
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FIG. 1: Evolution of distribution functions in free streaming cases, (a) Distribution functions at a 
cross section x = when Eq. ()lip is taken as an initial condition. The thick solid line represents the 
initial distribution function. The time advanced distribution function is given by the black dots, 
while the analytical solution is given by the dash-dotted curve, (b) Time advanced distribution 
function at x = vr (solid curve) and x = Sir (dash-dotted curve) when Eq. (ll3p is taken as an initial 
condition. 14 



1.0e-01 



1.0e-02 \ 




1 .Oe-06 



20 40 60 

t 

FIG. 2: Linear simulation results. Damping of the electric field strength \E\ is shown. Here, 
Vmax = 4:-0ve, and < X < Att. The mesh points are 8 and 32 for x and f , respectively. 
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FIG. 3: Nonlinear simulation results: (a) Time evolution of electric field strength \E\; (b) The 
initial distribution function (black solid curve) and the distribution function at t = 75 (red solid 
curve), both given at x = vr. 
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FIG. 4: An example of kappa distribution function: (a) Maxwellian and k = 1 function compared 
in regular scale; (b) k = 1 (red) , k = 2 (green), k = 5 (blue), and Maxwellian (black) compared 
in a logarithmic scale. 



17 




FIG. 5: Linear damping of electric field strength \E\ comparing the cases when Maxwellian (black) 
and K = 1 (red) and k = 2 distribution functions are taken as initial conditions. 
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FIG. 6: Effect of kappa function summarized: (a) Damping rate normalized by the absolute value 
of damping rate with Maxwellian; (b) real frequency normalized by the value of frequency with 
Maxwellian. 
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FIG. 7: The absolute value of d-uf{v) at the phase velocity. For the phase velocity u/k, the values 
of uj are employed from Fig. 6(b). The plotted values are normalized by that of the Maxwellian. 
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FIG. 8: Comparison of the numerical roots of the dispersion relation and the simulation results: 
(a) The damping rates 7 versus the wave vector k ;(b) the real frequencies w versus k. The roots 
from the (modified) dispersion relation are plotted as dash-dotted curves. The black, red, and 
green dash-dotted curves are for k = 2,3, and 6, respectively. Black circles are obtained from 
linear numerical simulation. Numerical simulation is done for k = 0.5,1.0, and 2.0 for k = 2,3, 
and 6. The simulation results are given at a fixed point x = vr. 
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FIG. 9: Electron distribution functions suggesting long time evolution up to t = 1500. The 
distribution functions are given at a fixed point x = n. (a) For a Maxwellian, (b) for a kappa 
distribution function (k = 2). The dash-dotted curves are for t = and the solid curves are for 
t = 1500. (c) Local expansion of (a) and (b) near the resonant phase velocities. Here, the black 
dash-dotted line at uj/k = 2.82 is for Maxwellian at t = 0, the black solid line at uj/k = 2.66 is for 
Maxwellian at t = 1500, the red dash-dotted line at oj/k = 2.32 is for kappa function at t = 0, and 
the red solid line at oj/k = 2.47 is for kappa function at t = 1500, respectively. The subscripts in 
the figure, max and kap stand for the Maxwellian and the kappa distribution function, respectively. 
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